#library(RMyCAL)
start_time=Sys.time()
library(GA)
source("function_modified_all.R")
resloc='./'

for (reti in seq(1)) {
  
Unsoda_code=read.table("./only_codes.csv",head=TRUE,sep=",")
nrow_code=nrow(Unsoda_code)
print(nrow_code)


#----------------------------------
#read all retention data
Unsoda_data_all=read.table("./lab_drying_h-t.csv",head=TRUE,sep=",")

#----------------------------------
#read all particle distribution data
Unsoda_data_particle_all=read.table("./particle_size.csv",head=TRUE,sep=",")


#--- for save data
Unsoda_data_retention=mat.or.vec(nrow_code,52*2+2)  # 790 samples, 52 max retention
Unsoda_fitted=mat.or.vec(nrow_code,41)  # for fitted parameters
index_nrow=mat.or.vec(nrow_code,2)

#--- make graphs
fname=paste("./fitted_Unsoda_test",".eps",sep="")   #!
postscript(file=fname, onefile=TRUE, paper="special",width=8.5, height=8.5, horizontal=F);
layout(matrix(seq(16), 4, 4, byrow = TRUE))
par(mai=c(0.3,0.3,0.3,0.3))

xl=c(0.0,0.8)
yl=c(0.01,80000)
suction_p=c(seq(0.01,0.09,by=0.01),seq(0.1,0.9,by=0.1),seq(1,300,by=2),seq(300,1200,by=50),seq(1200,80000,by=500))  # cm   pressure head

#---fit and make graph
# seq(nrow_code)
for (i in seq(nrow_code)) {
# for (i in seq(nrow_code)) {
  
  print(i)
  index=Unsoda_code[i,]
  #index=2362
  #DBconn=OpenDB_unsoda()  
  #Unsoda_data=Unsoda_sql(DBconn,index)
  #CloseDB(DBconn)
  Unsoda_data=Unsoda_data_all[Unsoda_data_all$code==index,]

  len=nrow(Unsoda_data)
  
  index_nrow[i,1]=index
  index_nrow[i,2]=len
  print(len)
  
  if ( len>=5 & sum(as.numeric(as.character(Unsoda_data$preshead))>=330)>=1 ){
    print("the index is:")
    print(index)
    print(Unsoda_data)
    
    #---particle size data
    #DBconn=OpenDB_unsoda()
    #Unsoda_data_particle=Unsoda_sql_particle(DBconn,index)
    #CloseDB(DBconn)
    Unsoda_data_particle=Unsoda_data_particle_all[Unsoda_data_particle_all$code==index,]

    #---save  Unsoda retention
    print(i)
    Unsoda_data_retention[i,1]=index
    Unsoda_data_retention[i,2]=len
    Unsoda_data_retention[i,seq(3,3+(len-1)*2,2)]=as.numeric(as.character(Unsoda_data$preshead))
    Unsoda_data_retention[i,seq(4,4+(len-1)*2,2)]=as.numeric(as.character(Unsoda_data$theta))
    
    Unsoda_fitted[i,1]=index
    #---------
    
    sample_h=as.numeric(as.character(Unsoda_data$preshead))
    sample_theta=as.numeric(as.character(Unsoda_data$theta))
    dat=data.frame(sample_number=as.numeric(len),h=as.numeric(sample_h),theta=as.numeric(sample_theta))
    
	
    #BC model, GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start BC model GA---")
      OLS_BC <- function(dat, thr, ths, psi_s, npar){
        attach(dat, warn.conflicts=F)
        
        theta_est=NULL
        for (i in seq(length(h)) ) {
          if (h[i]>psi_s) {
            #print(psi_s)
            #print(h[i])
            theta_est[i] =  thr+(ths-thr)* (psi_s/h[i])^npar
          } else {
            theta_est[i]=ths
          }
        }
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        #print(SSE)
        return(RMSE)
      }
      
      ga.bc <- ga(type='real-valued', min=c(0, 0.2, 1, 0.01), 
                  max=c(0.2, 0.8, 3000, 5), popSize=2, maxiter=2, #names=c('intercept', 'Wind', 'Temp'),
                  keepBest=T, fitness = function(b) -OLS_BC(dat, b[1],b[2],b[3],b[4]))#,  monitor=gaMonitor2)
      
      ga.model = summary(ga.bc)
      bc_solution=ga.bc@solution[1, ]
      bc_fitness=-ga.bc@fitnessValue[1]
      print(ga.model)
      bc_par=data.frame(theta_r=bc_solution[[1]],theta_s=bc_solution[[2]],psi_s=bc_solution[[3]],n=bc_solution[[4]])
      theta_bc=BC2(bc_par,suction_p)
      theta_h_bc=data.frame(theta=theta_bc,pressure=suction_p)
      #print(theta_h_bc)
      
      Unsoda_fitted[i,2:5]=bc_solution
      Unsoda_fitted[i,6]=bc_fitness
      
      print(bc_solution)
      print(bc_fitness)
    }
    
    #VG model, GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start VG model GA---")
      OLS_VG <- function(dat, thr, ths, alp, npar){
        attach(dat, warn.conflicts=F)
        theta_est = (thr+(ths-thr)/(1+(alp*abs(h))^npar)^(1-1/npar))
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        #print(SSE)
        return(RMSE)
      }
      
      ga.vg <- ga(type='real-valued', min=c(0, 0.2, 0.00033, 1.0), 
                  max=c(0.2, 0.8, 1, 6), popSize=2, maxiter=2, #names=c('intercept', 'Wind', 'Temp'),
                  keepBest=T, fitness = function(b) -OLS_VG(dat, b[1],b[2], b[3], b[4]) )#  , monitor=gaMonitor2)
      
      ga.model = summary(ga.vg)
      vg_solution=ga.vg@solution[1, ]
      vg_fitness=-ga.vg@fitnessValue[1]
      #print(ga.model)
      vg_par=data.frame(thr=vg_solution[[1]],ths=vg_solution[[2]],alp=vg_solution[[3]],n=vg_solution[[4]])
      theta_vg=vG(vg_par,suction_p)
      theta_h_vg=data.frame(theta=theta_vg,pressure=suction_p)
      
      Unsoda_fitted[i,7:10]=vg_solution
      Unsoda_fitted[i,11]=vg_fitness
      
      print(vg_solution)
      print(vg_fitness)
    }
    
    #KO(modified) model GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start modified KO GA---")
      OLS_ko <- function(dat, thr, ths , alp, n){
        attach(dat, warn.conflicts=F)
        theta_est = thr + (ths-thr)/2 - (ths-thr)/pi * atan((log(h)-alp)/n)
        if(sum(is.na(theta_est))>=1){
          theta_est[which(is.na(theta_est)==T)] <- ths
        }
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        return(RMSE)
      }
      
      ga.ko  <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                   upper=c(0.2, 0.8, 20, 20), popSize=2, maxiter=2, #names=c('intercept', 'Wind', 'Temp'),
                   keepBest=T, fitness = function(b) -OLS_ko(dat, b[1],b[2], b[3], b[4])  )#, monitor=gaMonitor2)
      
      #### summary of the ga with solution
      ga.model = summary(ga.ko)
      
      ko_solution=ga.ko@solution[1, ]
      ko_fitness=-ga.ko@fitnessValue
      
      ko_par=data.frame(thr=ko_solution[[1]],ths=ko_solution[[2]],alp=ko_solution[[3]],n=ko_solution[[4]])
      theta_ko=kom(ko_par,suction_p)
      theta_h_ko=data.frame(theta=theta_ko,pressure=suction_p)
      
      Unsoda_fitted[i,12:15]=ko_solution
      Unsoda_fitted[i,16]=ko_fitness
      
      print(ko_solution)
      print(ko_fitness) 
    }
    
    #lo(modified) model GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start modified HP GA---")
      OLS_lo <- function(dat, thr, ths , alp, n){
        attach(dat, warn.conflicts=F)
        theta_est = thr+(ths-thr)*0.5*(1-tanh((log(h)-alp)/n))
        if(sum(is.na(theta_est))>=1){
          theta_est[which(is.na(theta_est)==T)] <- ths
        }
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        return(RMSE)
      }
      
      ga.lo  <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                   upper=c(0.2, 0.8, 20, 20), popSize=2, maxiter=2, #names=c('intercept', 'Wind', 'Temp'),
                   keepBest=T, fitness = function(b) -OLS_lo(dat, b[1],b[2], b[3], b[4])  )#, monitor=gaMonitor2)
      
      #### summary of the ga with solution
      ga.model = summary(ga.lo)
      
      lo_solution=ga.lo@solution[1, ]
      lo_fitness=-ga.lo@fitnessValue
      
      lo_par=data.frame(thr=lo_solution[[1]],ths=lo_solution[[2]],alp=lo_solution[[3]],n=lo_solution[[4]])
      theta_lo=lom(lo_par,suction_p)
      theta_h_lo=data.frame(theta=theta_lo,pressure=suction_p)
      
      Unsoda_fitted[i,17:20]=lo_solution
      Unsoda_fitted[i,21]=lo_fitness
      
      print(lo_solution)
      print(lo_fitness) 
    }
    
	  #at(modified) model GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start modified AT GA---")
      OLS_at <- function(dat, thr, ths , alp, n){
        attach(dat, warn.conflicts=F)
        theta_est = thr + (ths-thr)/2 - (ths-thr)/pi * atan((log(h)-alp)/n)
        if(sum(is.na(theta_est))>=1){
          theta_est[which(is.na(theta_est)==T)] <- ths
        }
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        return(RMSE)
      }
      
      ga.at  <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                   upper=c(0.2, 0.8, 20, 20), popSize=2, maxiter=2, #names=c('intercept', 'Wind', 'Temp'),
                   keepBest=T, fitness = function(b) -OLS_at(dat, b[1],b[2], b[3], b[4])  )#, monitor=gaMonitor2)
      
      #### summary of the ga with solution
      ga.model = summary(ga.at)
      
      at_solution=ga.at@solution[1, ]
      at_fitness=-ga.at@fitnessValue
      
      at_par=data.frame(thr=at_solution[[1]],ths=at_solution[[2]],alp=at_solution[[3]],n=at_solution[[4]])
      theta_at=atm(at_par,suction_p)
      theta_h_at=data.frame(theta=theta_at,pressure=suction_p)
      
      Unsoda_fitted[i,22:25]=at_solution
      Unsoda_fitted[i,26]=at_fitness
      
      print(at_solution)
      print(at_fitness) 
    }
    
	  #gd model, GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start gd model GA---")
      OLS_gd <- function(dat, thr, ths, alp, n){
        attach(dat, warn.conflicts=F)
        gdx  <- function(t) atan(sinh(t))
        theta_est = thr + (ths - thr)/2 - (ths - thr)/pi *gdx((log(h)-alp)/n)
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        #print(SSE)
        return(RMSE)
      }
      
      ga.gd <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                  upper=c(0.2, 0.8, 20, 20), popSize=2, maxiter=2,
                  keepBest=T, fitness = function(b) -OLS_gd(dat, b[1],b[2], b[3], b[4]) )
      
      ga.model = summary(ga.gd)
      gd_solution=ga.gd@solution[1, ]
      gd_fitness=-ga.gd@fitnessValue[1]
      #print(ga.model)
      gd_par=data.frame(thr=gd_solution[[1]],ths=gd_solution[[2]],alp=gd_solution[[3]],n=gd_solution[[4]])
      theta_gd=gdm(gd_par,suction_p)
      theta_h_gd=data.frame(theta=theta_gd,pressure=suction_p)
      
      Unsoda_fitted[i,27:30]=gd_solution
      Unsoda_fitted[i,31]=gd_fitness
      
      print(gd_solution)
      print(gd_fitness)
    }
    
    #e2 model GA algorithm
    #pjli 2022.03.19.change constant coef
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start e2 GA---")
      OLS_e2 <- function(dat, thr, ths, alp, n){
        attach(dat, warn.conflicts=F)
        theta_est = ths -(ths - thr)*(1 + exp(-(log(h)-alp)/n))^(-1)
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        return(RMSE)
      }
      
      ga.e2 <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                  upper=c(0.2, 0.8, 20, 20), popSize=1000, maxiter=700, 
                  keepBest=T, fitness = function(b) -OLS_e2(dat, b[1],b[2], b[3], b[4])  )#, monitor=gaMonitor2)
      
      #### summary of the ga with solution
      ga.model = summary(ga.e2)
      #print(ga.model)
      e2_solution=ga.e2@solution[1, ]
      e2_fitness=-ga.e2@fitnessValue  
      
      e2_par=data.frame(thr=e2_solution[[1]],ths=e2_solution[[2]],alp=e2_solution[[3]],n=e2_solution[[4]])
      theta_e2=e2m(e2_par,suction_p)
      theta_h_e2=data.frame(theta=theta_e2,pressure=suction_p)
      
      Unsoda_fitted[i,32:35]=e2_solution
      Unsoda_fitted[i,36]=e2_fitness
      
      print(e2_solution)
      print(e2_fitness)  
    }
    
    #sq model GA algorithm
    #-----------
    if(T) {
      print(paste('index=',index,'     num=',i,'     tempo=',format(i*100/nrow_code, digits = 3)))
      print("---start sq GA---")
      OLS_sq <- function(dat, thr, ths , alp, n){
        attach(dat, warn.conflicts=F)
        theta_est = thr + (ths-thr)/2 - (ths-thr)/2 * (log(h)-alp)/n * (1 + ((log(h)-alp)/n)^2)^(-1/2)
        if(sum(is.na(theta_est))>=1){
          theta_est[which(is.na(theta_est)==T)] <- ths
        }
        SSE= t(theta-theta_est) %*% (theta-theta_est)
        df=length(theta)
        RMSE=sqrt(SSE/df)
        detach(dat)
        return(RMSE)
      }
      
      ga.sq  <- ga(type='real-valued', lower=c(0, 0.2, 0, 0.01), 
                   upper=c(0.2, 0.8, 20, 20), popSize=2, maxiter=2, 
                   keepBest=T, fitness = function(b) -OLS_sq(dat, b[1],b[2], b[3], b[4])  )
      
      #### summary of the ga with solution
      ga.model = summary(ga.sq)
      
      sq_solution=ga.sq@solution[1, ]
      sq_fitness=-ga.sq@fitnessValue
      
      sq_par=data.frame(thr=sq_solution[[1]],ths=sq_solution[[2]],alp=sq_solution[[3]],n=sq_solution[[4]])
      theta_sq=sqm(sq_par,suction_p)
      theta_h_sq=data.frame(theta=theta_sq,pressure=suction_p)
      
      Unsoda_fitted[i,37:40]=sq_solution
      Unsoda_fitted[i,41]=sq_fitness
      
      print(sq_solution)
      print(sq_fitness) 
    }
    
    
    print(head(Unsoda_data_retention))
    
    colnames(Unsoda_fitted)=c("code","thr_BC","ths_BC","psi_BC","n_BC","rmse_BC",
    								 "thr_VG","ths_VG","alp_VG","n_VG","rmse_VG",
    								 "thr_KO","ths_KO","alp_KO","n_KO","rmse_KO",
    								 "thr_HP","ths_HP","alp_HP","n_HP","rmse_HP",
    								 "thr_AT","ths_AT","alp_AT","n_AT","rmse_AT",
    								 "thr_GD","ths_GD","alp_GD","n_GD","rmse_GD",
    								 "thr_SG","ths_SG","alp_SG","n_SG","rmse_SG",
    								 "thr_CA","ths_CA","alp_CA","n_CA","rmse_CA")
    name1=paste(resloc,"/Unsoda_fitted",reti,".csv",sep="")
    write.table(Unsoda_fitted,name1,row.name=FALSE, quote=FALSE,col.name=TRUE,sep=',')
    name2=paste(resloc,"/Unsoda_retention_data",reti,".csv",sep="")
    write.table(Unsoda_data_retention,name2,row.name=FALSE, quote=FALSE,col.name=FALSE,sep=',')
  }
}

end_time=Sys.time()
print(paste('start time',start_time))
print(paste('end time  ',end_time))
}
